clear all; clear global; clc;

Cal.start_time = datestr(now, 'HH:MM:SS');

load us_rat_smooth

Cal.init_year = 1975;
Cal.mid_year  = 1981;

us_rat_1975 = 100*us_rat_1975/sum(us_rat_1975);
us_rat_1981 = 100*us_rat_1981/sum(us_rat_1981);

Cal.y_init = 100/33*ones(33,1)/100;
Cal.y_init = us_rat_1975/100;

Cal.mu_init = us_rat_1975(1:32)/100;
Cal.mu_last = us_rat_1981(1:32)/100;

Cal.RnD_int_A = 1.75/100;
Cal.RnD_int_B = 1.96/100;

Cal.Growth_A = 0.55/100;
Cal.Growth_B = 1.73/100; 

Cal.Entry_A = 0.10;
Cal.EXtoY_A = 0.07;

Cal.Prft_rat = .06;

Cal.tau_A = 0.053;
Cal.tau_B = 0.038;

Cal.L_A   = 1;
Cal.L_B   = 1;

% ====== Bounds for each Parameter to Calibrate =========

alfa_Low   = 0.001;             
alfae_Low  = 0.001;            
phi_Low    = 0.001;            
dlt_Low    = 0.001;    
lmd_Low    = 0.0001;           
kapa_Low   = 0.0001;           

alfa_High   = 5;             
alfae_High  = 5;            
phi_High    = 100.999;       
dlt_High    = 0.5; 
lmd_High    = 0.25;          
kapa_High   = 0.5;          

Cal.Low = [alfa_Low alfa_Low alfae_Low alfae_Low phi_Low dlt_Low lmd_Low kapa_Low];
Cal.High= [alfa_High alfa_High alfae_High alfae_High phi_High dlt_High lmd_High kapa_High];

Cal.Low = [alfa_Low alfa_Low alfae_Low alfae_Low phi_Low lmd_Low kapa_Low];
Cal.High= [alfa_High alfa_High alfae_High alfae_High phi_High lmd_High kapa_High];

%% ============= Guess ===============

load diff_new 
guess_new = Param0;

Eq_new = @(M)Equilibrium(M,Cal);

opt_min = optimset('MaxFunEvals',20000,'MaxIter',20000,'TolFun',.01,'MaxTime',3600);

[diff_min,f_min,exitflag,output] = fminsearch(Eq_new,guess_new,opt_min);

Param = exp(diff_min)./(1+exp(diff_min)).*Cal.High+(1-exp(diff_min)./(1+exp(diff_min))).*Cal.Low;
save Param_cal Param
save diff_cal diff_min